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We investigate the H 3 S phase of sulphur hydride under high pressure ~ 200 GPa by means of 
ab-initio calculations within the framework of the density-functional theory (DFT) with the PBEO 
hybrid exchange-correlation (i?xc) approximation. The choice of Exc has the largest effect on the 
calculated electron-phonon coupling (EPC) matrix elements; the high pressure equation of state 
and phonon frequencies are only slightly modified. Mode-dependent EPC correction factors are 
determined from PBEO using a frozen-phonon supercell approach, while standard density-functional 
perturbation theory is used to determine the EPC with PBE generalized-gradient approximation 
Exc- Our principle finding is that the calculated PBEO Tc is enhanced by 25% compared to PBE. This 
is similar in magnitude, but in opposite direction, to the proposed suppression of Tc by anharmonic 
effects [Errea et al, Phys. Rev. Lett. 114, 157004 (2015)]. Our calculations demonstrate the 
importance of considering exchange-correlation approximations for calculations of superconducting 
properties for this class of materials. 
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The quest for the holy grail of high-pressure physics, 
metallic hydrogen, has continued to attract the interest 
of experimentalists and theorists since Ashcroft proposed 
that the new phase should exhibit superconductivity with 
Tc ~ 270 K.ffl A recent focus has been on the hydrides 
where reduced metallization pressures are expected.^ Sul¬ 
phur hydrides have attracted a great deal of attention 
due to theoretical predictiond^^^ of high ~ 200 K un¬ 
der high pressure. This was supported by two experi¬ 
mental reports!^ superconductivity was first attributed 
to H 2 S in Ref. 0 but Drozdov et alW re-analyzed their 
measurements, which lead to « 203 K, now ascribed 
to the H 3 S phase. Theoretically, Li et proposed a 
metallic phase of hydrogen sulfide H 2 S, potentially su¬ 
perconducting with a maximum Tc ^ 80 K when sub¬ 
jected to a pressure of 160 GPa. The transition temper¬ 
ature was estimated by applying the Alien-Dynes modi¬ 
fied McMillan equatiorPSl and assuming electron-phonon 
comling as the source of superconductivity. Duan et 
al^ performed a systematic investigation of the pressure- 
dependent (H 2 S) 2 H 2 phase diagram and concluded that 
the cubic Im3m H 3 S structure was the most stable phase 
at pressures above 180 GPa, with Tc ~ 200 K at 200 
GPa. Papaconstantopoulos et al^ applied the Gaspari- 
Gyorffy theorjff^ and argued that superconductivity in 
H 3 S arose mostly from the coupling between the electrons 
and the H vibrations, whereas the role of the sulphur is to 
stabilize the hydride at high pressures via hybridization. 
A similar hypothesis was put forward by Bernstein et alW 
who proposed that the transport mechanism in the high- 
pressure H 3 S was the same as in MgB 2 . A substantially 
higher T^, in the hydride could be explained with the con¬ 
siderably smaller atomic masses of the constituents. It 
has also been argued that the features of the calculated 
electron-phonon spectrum in H 3 S are near optimum for 


a high Tc.^ Errea et claimed, however, that the har¬ 
monic approximation, used in the previous calculations, 
overestimated the electron-phonon coupling and conse¬ 
quently the predicted Tc. Other related candidate ma¬ 
terials, such as SeH 3 , have also been investigated.^i^ By 
contrast, Hirsch and MarsigliJi^ have argued against the 
conventional electron-phonon mechanism and in favor of 
an electron-hole mediated mechanism. 

The previous theoretical investigation^^HZl were based 
on ab-initio calculations that used the semilocal PBE!^ 
generalized-gradient approximation (GGA). Standard lo¬ 
cal or semilocal DFT approaches sometimes fail to pre¬ 
dict the correct Tc in conventional electron-phonon su¬ 
perconductors because of an insufficient treatment of the 
exchange-correlation effects.^i^lti^ In this paper we ad¬ 
dress the effect of the choice of DFT exchange-correlation 
functional on the predicted EPC strength. Hybrid DFT 
(HDFT) methods often improve the p redict ed proper¬ 
ties of materials with hydrogen bonds.^^SHI] Jn HDFT 
a fraction of the Hartree-Fock exact-exchange term is 
added to E^c^^ a result, HDFT calculations are typ¬ 
ically at least an order of magnitude computationally 
more demanding than local or semilocal Fxc- To cal¬ 
culate the effective Kohn-Sham potential for HDFT, it 
is necessary to evaluate A^k x Nocc six-dimensional spa¬ 
tial integrals, where A^k and Nqcc are the number of k- 
points and occupied states, respectively. By contrast, 
only three-dimensional integrals are required for local or 
semilocal Axe- For this reason, the extension and appli¬ 
cation of DFT perturbation theory (DFTPT) to calcu¬ 
late electron-phonon coupling coefficient^^ becomes very 
challenging in HDFT. Instead, we adapted the approxi¬ 
mation scheme introduced by Lazzeri et and Yin et 
for GW and HDFT, as described below. 

The mode-dependent EPC strengths Aqi, are given by 
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where the phonon frequencies Wq^ and eigenvectors Uqj^ 
are obtained by diagonalizing the phonon dynamical ma¬ 
trix; ek,i, ^k,i are the single-electron energies and eigen¬ 
functions, Niep) is the density of states at the Fermi en¬ 
ergy ef, and Mj is the mass of the /th atom. The EPC 
matrix element ('0k,i| I0k-i-q,j) depends on the 

change of the self-consistent Kohn-Sham potential with 
respect to the phonon displacement. All of the quanti- 
ties in Eq. Q can be calculated using DETPlt^SHUl for 
local or semi-local E'xc, and this functionality is available 
in standard codes such a Quantum Espresso.!^ Unfor¬ 
tunately, a combination of HDET and DEPT would be 
computationally too demanding and has not yet been im¬ 
plemented for solids, to our k nowledge. Instead, we use 
an alternative approach. b^P^ I Supercells and the frozen- 
phonon method are used to obtain the phonon dynam¬ 
ical matrix and the resulting pairs of Wqi,, Uq,^ on a q- 
grid commensurate with the supercell.l^ Small but finite 
phonon displacements can then be introduced to obtain 
phonon-induced Kohn Sham potentials, which can be ex¬ 
pressed to first order as: U{r} + (dV{R}/dR/) • The 
change Ae of a single-electron energy is approximated as: 

(^) 


The difference Ae describes band splitting due to the 
lifted degeneracy in the presence of the phonon. For the 
application of the approximation for the calculation 
of Aqy only the splittings nearest to the Fermi level are 
relevant. An example is presented in Fig. for ImSm 
H 3 S in the presence of phonon mode = 6 at q = 0 , 
calculated with PBE E^c- 

The band splittings are determined for both DFT and 
HDFT, and respectively. Using Eq. (1^ 

in Eq. Q for DFT and for HDFT yields the following 
approximation to estimate the HDFT electron-phonon 
coupling coefficients: 


\ HDFT ^ \DFT p 
^qjy ~ -^qv Jqi^J 
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where the correction factor is given by: 


_^DFT(,DFT) /^DFXy /^.HDFXy 
_^HDFT (gHDFT) (^ccHDFT j ^ ^,DFT j ' 


This procedure could be repeated, in principle, to de¬ 
termine for each phonon mode throughout the 

Brillouin zone (BZ), using appropriate q-commensurate 
supercells. We note, however, that the Alien-Dynes equa- 
tion for the superconducting transition temperaturelSl 


Tc = 


CJiog 

TJ 
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/ -1.04(1-f A) A 
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FIG. 1. Band splitting due the presence of a phonon. PBE 
band structures of ImSm H 3 S are shown: solid lines corre¬ 
spond to the undistorted crystal, and dashed lines are for 
the atoms displaced according to the phonon mode 1 / = 6 at 
q = 0. The inset presents the details of the band splitting 
close to the Fermi energy near the F point due to the lifted 
degeneracy caused by the presence of the phonon. 


depends on the total electron-phonon coupling coefficient 
A, which is given by an integral over the BZ and sum over 
phonon branches: A = (1 /Dbz) 1 where Obz 

is the volume of the BZ; fi* and wiog are discussed below. 
Therefore we determine the center of gravity q* in the 
irreducible wedge of the Brillouin zone (IBZ) with respect 
to the AqJ"’’ for each DFT phonon branch f: 


/jBzd^^ADFT ■ 


( 6 ) 


For each branch we calculate the correction factor (4) 
/q-y and use this single value for all q-points within a 
given branch v: 


X HDFT ^ \DFT p 
^qjy ~ '^qu Jq^v 


( 7 ) 


The calculations were carried out by applying the 
Quantum EspresscP^ code. The DFT exchange- 
correlation po tentia l was PBE, whereas for HDET 
we used PBEO.^^HIIl xhe bare electron-ion interactions 
were described with t he no rm-conserving Goedecker- 
Hartwigsen-Hutter-Teteil22HU pseudopotentials. The 
planewave and charge-density cut-off parameters were 
set to 476 eV and 1904 eV, respectively. A Monkhorst- 
PaclP^ 16 X 16 X 16 k-point grid was used for the BZ 
integration primitive cell calculations. PBEO phonon 
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FIG. 2. The calculated total energy with respect to the lattice 
parameter a for PBE (+) and PBEO (x) fitted with the Birch- 
Murnaghan equation of state, from which the pressure was 
calculated. The solid lines are for PBE, the dashed lines are 
for PBEO. 

spectra were calculated using the frozen-phonon method 
as implemented in the Phonopji^ code, whereas the 
DFT electron-phonon coefficients were determined using 
DFPT, which is a part of the Quantum Espresso pack¬ 
age. The total electron-phonon coupling coefficient 
was obtained by using a 4 x 4 x 4 q-point grid in the BZ. 
The force constants for the frozen-phonon method were 
obtained from a 2 x 2 x 2-supercell calculation. 

Fig. 2 presents the calculated total energies as a func¬ 
tion of the lattice parameter a fitted with the Birch- 
Murnaghan equation of state.l^ The difference between 
PBE and PBEO is most pronounced at zero pressure, 
reflecting the ~ 2% difference in equilibrium lattice pa¬ 
rameters: 3.66 A and 3.58 A, respectively. Near 200GPa, 
however, PBE and PBEO yield ve^ similar lattice param¬ 
eters, which match the published^ value 2.984A. We used 
this value for the subsequent calculations of the phonon- 
related properties. 

The effect of HDFT is larger for phonon frequencies 
near 200 GPa, as shown in Fig. 3. The largest differences 
are for the optic modes, especially near the N point, while 
differences for the acoustic modes are small. PBE EPG 
strengths throughout the Brillouin zone are depicted in 
Fig. 4. In contrast to the case of BaBiOa from Ref. [HI 
there are no distinguished points with particularly large 
values of This provides some justification for our 

introduction here of the center of gravity point q* in 
Eq. (5). Using this approximation, the number of HDFT 
fqi, evaluations is reduced from 96 (the number of q- 
points in the irreducible BZ times the number of bands) 
to only 12 (the number of bands) in Eqs. (4) and (6). For 
all phonon branches, the center of gravity point is well- 
approximated by either q* = F or q* = (0,0,1/2) 27r/a. 
The band splittings Ae can therefore be obtained from 



FIG. 3. Galculated phonon dispersion for the PBE (solid 
lines) and PBEO (dashed lines) exchange-correlation poten¬ 
tials. 



EIG. 4. Electron-phonon q-dependent PBE coupling 
strengths: the size of the circles, superimposed on the phonon 
dispersion curves, is proportional to the electron-phonon cou¬ 
pling coefficient AqJ^. 


primitive cell or 1 x 1 x 2-supercell calculations, respec¬ 
tively. 

A suitable indicator for the overall influence 
of the exchange-correlation effects on the electron- 
phonon coupling is the Eliashberg spectral function 
a^F{w) oc J(l/fiBz)\,^cOq^S(uj - ojq^)d^q which is 
plotted in Fig. 5 (using gaussians with width tr = O.STHz 
to represent the ^-functions). The difference in peak po¬ 
sitions reflects the changes in phonon dispersion in Fig. 3. 
The highest peaks are located in the middle of the fre¬ 
quency range, which is in agreement with the distribution 
of the Aq®^ magnitudes in Fig. 4. The PBEO peaks for 
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FIG. 5. The calculated Eliashberg spectral function a^F{uj), 
using PBE (solid lines) and the PBEO (dashed lines) 
exchange-correlation potentials. 



A 

taiog [K] 

TclK] 

PBE 

1.76 

1657 

217-201 

PBEO 

2.18 

1773 

270-253 


TABLE I. A comparison between the PBE and PBEO calcu¬ 
lated quantities, which appear in the Alien-Dynes equation 
1^, and between the resulting Tc for the retarded Coulomb 
repulsion jj.* = 0.1 — 0.13. 


frequencies up to 50 THz are higher than the correspond¬ 
ing PBE peaks, which is due to the enhanced phonon- 
induced band splitting in PBEO [Eq. (4)]. The situation 
is the opposite at the frequencies higher than 50 THz, 
where the PBE peaks are higher. 

The quantitative results are summarized in Ta¬ 
ble I. The superconducting transition temperature 
Tc in the Alien-Dynes equation © depends not 
only on the electron-phonon coupling coefficient 
A, but also on the logarithmic average frequency 
<^iog = (l/-^^Bz) / \i^^og{ujq^)d^q and the retarded 

Coulomb repulsion /i*. Here, we simply adopt literature 
values for p,*.Iii^From visual inspection of Fig. 3 it is not 
evident whether PBEO yields higher or lower phonon 
frequencies on average than PBE. Indeed, the calculated 
wiog values differ by only 7%. Thus, the change in 
phonon frequencies has only a small effect on the PBEO 
enhancement of A by 24%. A comparable enhancement 
of A was also found for graphene/graphit^i^ and Cgo 
molecule.lA^ As shown in the Table, the calculated 
PBEO Tc, using the Alien-Dynes equation, is higher 
by about 25% compared to PBE. The influence of the 
exchange-correlation effects on the /r* is complex, and 
it is not in the focus of the present paper. DFT and 


Hartree-Fock (and thus also HDFT) wave functions dif¬ 
fer only slightly,!^ so the double Fermi-surface averaged 
electron-electron Coulomb interactions will be similar. 
Since these determine fj. and it is reasonable to 

compare predictions of DFT and HDFT using the same 

. On the basis of the presented results it is clear that 
the predicted Tc would be enhanced for any reasonable 
change in fF after switching from PBE to PBEO. 

The context for our results is the theory of conven¬ 
tional electron-phonon superconductivity as described by 
DFT. In H 3 S and related sulfides, the normal-state Fermi 
liquid quasiparticles have been described by DFT band 
structures and wave functions, using PBE Txc. The 
DFT approximation is often quite good and is widely 
used, although there are no guarantees regarding the 
accuracy of the related one-particle Kohn-Sham eigen¬ 
states. Phonon excitations are commonly described in 
the harmonic approximation, using DEPT or DFT su¬ 
percell calculations. These are ground state properties, 
within the Born-Oppenheimer approximation, so DFT is 
well justified. Using DFT for the no n-adia batic EPC is 
also justified within Migdal’s theoremj ^^ l ^^ l and these are 
treated, as described above. The superconducting tran¬ 
sition temperature Tc is then usually estimated using the 
McMillan or Alien-Dynes equations, although t he super- 
conducting anisotropic gap equation can be usecpH^ESEZl 
by means of SCDFT. The SCDFT calculation with har¬ 
monic phonon^i^ yielded Tc =180 K for H 3 S, in agree¬ 
ment with experiment. While SCDFT provides a sound 
footing for the DFT treatment of superconductivity, it 
does not in itself improve the electronic band structure 
and phonon properties. The effects of anharmonicity 
were studied by Errea et alW They found that anhar¬ 
monicity suppressed A by 30%; T^ was suppressed by 34 
and 56 K, respectively, using the McMillan equation or 
the isotropic Migdal-Eliashberg equations. All of these 
studies were based on PBE E^c- Here we have shown 
that the choice of exchange-correlation potential can have 
effects of similar magnitude on the EPC and Tc- The 
present PBEO calculations modify the electronic quasi¬ 
particles in the theory, the phonon excitations, and the 
EPC. For Tc, we relied on the Alien-Dynes equations with 
our calculated PBEO A and phonon frequencies. 

In summary, we examined the influence of exchange- 
correlation effects on the electron-phonon coupling in 
cubic /m3m H 3 S under high pressure. We introduced 
the electron-phonon coupling center-of-gravity point q*, 
modifying the approach used in Ref. HU We calculated 
a 25% enhancement of Tc for HDFT compared to PBE 
predictions. This is similar in magnitude, but in oppo¬ 
site direction, to the proposed suppression of Tc by an- 
harmonic effects in Ref. [U Our results demonstrate the 
importance of considering exchange-correlation approx¬ 
imations for calculations of superconducting properties 
for this class of materials. 
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